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ABSTRACT 

Crosshole  radar  tomography  is  increasingly  being  used  to  characterize  the  shallow  subsur¬ 
face  and  to  monitor  hydrologic  processes.  Although  tomographic  inversion  provides  a  subsurface 
model,  confidently  interpreting  the  resulting  image  can  be  challenging.  We  conducted  a  simple 
modeling  study  to  better  understand  the  capabilities  and  limitations  of  tomographic  inversion.  We 
start  with  a  known  earth  model,  simulate  a  tomography  experiment,  and  invert  the  synthetic  data. 
We  investigate  the  effects  of  straight  and  curved  ray  approximations  to  wave  propagation,  regular¬ 
ization,  grid  size,  and  starting  model.  We  also  investigate  the  effects  of  limited  ray  coverage 
through  the  earth  model  and  noise  in  the  data.  Understanding  the  effects  of  these  different  meth¬ 
ods  and  parameterizations  will  help  us  place  confidence  limits  on  modeled  features  to  more  accu¬ 
rately  reflect  our  knowledge  of  the  subsurface. 

INTRODUCTION 

Crosshole  radar  tomography  is  increasingly  being  used  to  characterize  the  shallow  subsur¬ 
face  and  to  monitor  hydrologic  processes  such  as  infiltration  through  the  vadose  zone  and  solute 
transport  through  the  saturated  zone.  While  collecting,  processing,  and  inverting  radar  tomogra¬ 
phy  data  is  not  difficult,  interpreting  the  resulting  images  with  confidence  can  be  a  challenge.  Are 
parameter  estimates  accurate,  or  are  they  artifacts  of  the  inversion  process?  What  constraints 
should  be  used  in  the  inversion  process  to  improve  results?  Methods  to  determine  the  quantitative 
resolution,  uncertainty,  and  reliability  of  tomographic  images  are  greatly  needed. 

As  a  first  step  towards  addressing  these  issues,  we  conducted  a  simple  modeling  study.  We 
start  with  a  known  earth  model,  simulate  a  tomography  experiment,  and  invert  the  synthetic  data. 
By  comparing  the  inversion  results  with  the  starting  model,  we  can  improve  our  understanding  of 
the  capabilities  and  limitations  of  different  tomographic  inversion  methods. 

We  used  two  different  inversion  algorithms:  one  straight  ray  and  one  curved  ray.  We  also 
investigated  the  effects  of  including  different  types  of  constraints  in  the  inversion,  as  well  as  the 
effects  of  ray  coverage,  uncertainty  in  traveltime  picks,  model  discretization,  and  starting  model. 
Good  references  on  tomography  and  inverse  theory  include  Menke  (1989)  and  Snieder  and 
Trampert  (1999). 

The  earth  model  that  we  selected  is  representative  of  the  types  of  features  and  electrical 
property  variations  that  we  expect  at  the  Boise  Hydrogeophysical  Research  Site  (BHRS),  the 
focus  of  an  intense  investigation  of  geophysical  and  hydrologic  parameter  distribution  in  a  hetero¬ 
geneous,  coarse,  alluvial  aquifer  (Barrash  et  al.,  1999;  Clement  et  al.,  1999;  Peterson  et  al.,  1999). 
The  synthetic  modeling  results  in  this  paper  provide  insight  into  the  capabilities  and  limitations  of 
different  tomography  methods  and  the  level  of  detail  that  can  be  reliably  interpreted  from  tomo¬ 
graphic  images  of  the  site. 
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BACKGROUND  THEORY 

Matrix  algebra  provides  a  convenient  way  to  express  forward  and  inverse  modeling  prob¬ 
lems  such  as  travel  time  tomography.  The  forward  problem  is: 

Gm  =  d, 

where  d  is  the  data  (vector),  G  is  the  kernel  function  (matrix),  and  m  is  the  model  (vector).  The 
data  are  the  projection  of  the  model  by  the  kernel  function.  The  kernel  function  represents  the 
physics  of  the  problem,  including  boundary  conditions  and  differential  equations. 

A  solution  to  the  inverse  problem  is: 

m  =  G'1d. 

The  inversion  process  involves  computing  the  inverse  of  matrix  G  and  then  multiplying  this 

matrix  by  the  data  to  compute  the  model.  Unfortunately,  G'1  can  be  difficult  to  compute,  particu¬ 
larly  if  the  matrix  is  ill-conditioned  (small  data  errors  cause  large  model  changes),  ill-posed  (mix- 
determined)  (Snieder  and  Trampert,  1999),  or  large  (too  many  parameters  for  available  computer 
memory).  In  these  cases,  a  weighted,  damped,  least  squares  approach  is  often  used  to  find  a  solu¬ 
tion  (Menke,  1989): 

mest  =  <m>  +[GTWeG  +  X2  Wm]’1GTWe[d-G<rn>] 

where  mest  is  the  best  fitting  model,  <m>  is  the  starting  model,  d  and  G  are  the  data  and  the  kernel 
as  before,  Wm  is  the  regularization  or  model  weighting  matrix,  We  is  the  data  weighting  matrix, 

and  X  is  the  weighting  factor  between  prediction  error  (overdetermined)  and  solution  length 
(underdetermined).  The  parameter  X  determines  how  much  the  data  influences  the  model  versus 
how  much  the  model  is  constrained  by  the  regularization.  For  X  =  0,  the  solution  depends  only  on 
the  data.  For  large  X  values,  the  solution  depends  more  on  the  regularization. 

The  least  squares  approach  minimizes  the  L2  norm  to  determine  the  optimal  solution,  i.  e., 

where  e,  are  the  differences  between  the  observed  and  calculated  data.  In  the  tomography  prob¬ 
lem,  the  data  are  first  arrival  times  and  we  try  to  minimize  the  difference  between  the  calculated 
and  observed  times. 

Tomographic  data  sets  often  consist  of  many  measurements.  The  result  is  that  matrices  can 
be  large,  sparse,  and  difficult  to  invert  directly.  A  variety  of  computational  methods  have  been 
developed  to  implement  matrix  inversions,  including  ART  (e.  g.,  Peterson  et  al,  1985),  SIRT  (e. 
g.,  Tweeton,  1988),  and  LSQR  (Paige  and  Saunders,  1982),  which  are  all  iterative  solvers.  We  use 
a  straight  ray  tomographic  inversion  using  SIRT  (Tweeton,  1988)  and  a  curved  ray  inversion  using 
LSQR  (Aldridge  and  Oldenburg,  1993). 

FORWARD  MODELING 

To  test  different  crosshole  traveltime  tomography  algorithms,  we  used  a  finite-difference 
approximation  to  the  eikonal  equation  (Vidale,  1990)  to  simulate  a  crosshole  radar  tomography 
experiment.  The  eikonal  equation  calculates  the  first  arrival  travel  time  of  a  propagating  wave 
through  a  gridded  velocity  field.  Figure  1  shows  the  input  velocity  model.  The  model  is  5  m  by  10 
m  with  a  horizontal  and  vertical  grid  spacing  of  0.05  m.  We  included  several  layers,  some  with  lat¬ 
eral  velocity  changes,  to  simulate  features  that  we  believe  exist  at  the  BHRS.  The  input  velocity 
model  consists  of  high  and  low  velocity  zones  and  thick  and  thin  layers  to  test  the  inversion 
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schemes.  The  upper,  layer,  2  m  thick,  represents  the  vadose  zone  with  a  velocity  of  0.140  m/ns. 
Velocities  in  the  saturated  zone  (2  -  10  m  depth)  range  from  0.070  to  0.095  m/ns.  Layers  are  hori¬ 
zontal  and  range  in  thickness  from  0.5  to  2.0  m.  The  layer  at  4.0  to  5.5  m  depth  includes  a  low 
velocity  inclusion  in  the  center  of  the  model,  between  2.0  and  3.5  m  on  the  horizontal  scale;  at  the 
edges  of  the  model,  the  velocity  is  0.090  m/ns  while  in  the  center  the  velocity  is  0.070  m/ns. 
Another  lateral  velocity  change  is  included  in  the  layer  between  7.0  and  8.0  m  depth.  Close  to  the 
wells,  the  velocity  is  0.080  m/ns  while  in  the  center,  between  2.0  and  3.5  m  on  the  horizontal 
scale,  the  velocity  is  0.090  m/ns.  These  horizontal  velocity  changes  are  included  to  determine  the 
ability  of  the  inversion  methods  to  image  heterogeneity  in  the  subsurface. 

We  simulated  a  traveltime  tomography  experiment  in  two  wells  spaced  3.5  m  apart.  Travel 
times  were  computed  between  40  shot  locations  (0.25  m  vertical  spacing)  and  41  receiver  loca¬ 
tions  (also  0.25  m  vertical  spacing). 

We  used  two  forward  modeling  algorithms  to  invert  the  synthetic  crosshole  traveltime 
data.  The  first  algorithm  uses  the  same  forward  model  scheme  as  was  used  to  generate  the  syn¬ 
thetic  data  (Aldridge  and  Oldenburg,  1993).  This  algorithm  can  simulate  ray-bending  at  velocity 
contrasts.  We  also  used  a  straight-ray  forward  modeling  algorithm  for  comparison  (Tweeton, 
1988).  In  both  cases,  the  simulated  tomographic  data  consists  of  1640  first  arrival  times. 

RESULTS  OF  INVERSION 
Straight  Ray  versus  Curved  Ray  Simulators 

Figure  1  shows  the  input  model  and  a  comparison  of  the  straight  ray  and  curved  ray  inver¬ 
sions  using  noise-free  data.  Velocity  values  were  determined  on  a  0.1  m  grid,  and  then  upscaled  to 
0.25  m  for  display.  The  straight  ray  inversion  smooths  the  velocity  in  each  pixel  using  an  average 
of  its  near  neighbors  to  compensate  for  the  fine  grid  spacing.  The  curved  ray  method  uses  regular¬ 
ization  to  condition  the  problem;  in  this  simulation,  we  used  a  flatness  constraint  with  10: 1  hori¬ 
zontal  to  vertical  regularization. 

Results  show  that  the  curved  ray  simulator  is  better  able  to  locate  boundaries  than  the 
straight  ray  simulator.  The  curved  ray  model  has  a  sharp  velocity  change  from  about  0.140  m/ns 
in  the  upper  2  m  to  velocities  less  than  0.100  m/ns  in  the  lower  8  m  of  the  model.  This  boundary  is 
smeared  in  the  straight  ray  model.  Both  models  also  provide  a  hint  of  the  alternating  high  and  low 
velocities  with  depth.  Again,  the  boundaries  are  more  distinct  in  the  curved  ray  model.  The  mod¬ 
els  also  show  lateral  velocity  changes  in  the  center  of  the  domain,  although  the  boundaries  of  the 
inclusions  are  not  particularly  well  resolved. 

How  do  the  two  inversion  solutions  compare  to  the  true  velocity  field?  The  straight  ray 
model  has  a  velocity  range  of  0.059  to  0.150  m/ns,  somewhat  greater  than  the  true  range,  with  an 
average  misfit  error  of  0.005  ns  and  an  RMS  misfit  error  of  0.102  ns.  The  range  of  error  values  for 
the  straight  ray  model  is  -0.509  to  0.539  ns.  The  curved  ray  model  has  a  velocity  range  of  0.069  to 
0.157  m/ns,  an  average  misfit  error  of  -0.012  ns,  and  an  RMS  error  of  0.239  ns.  The  misfit  range 
for  the  curved  ray  model  is  -0.863  to  1.836  ns.  Both  methods  have  an  average  error  near  zero,  but 
the  curved  ray  method  has  a  larger  RMS  error  indicating  a  larger  variance  in  the  misfit  error. 

To  gain  insight. into  how  the  two  inversion  methods  partition  error  along  ray  paths,  we 
generated  plot  of  ray  density  for  both  methods.  Ray  density  is  defined  as  the  total  raypath  length 
per  slowness  cell  divided  by  the  grid  interval.  Figure  2  show  the  ray  density  of  the  two  inversion 
methods.  For  the  straight  ray  inversion,  more  rays  pass  through  the  center  of  the  model  than  at  the 
edges.  The  ray  coverage  decreases  systematically  away  from  the  model’s  center.  This  coverage 
difference  will  cause  the  inversion  algorithm  to  assign  errors  to  the  less  densely  sampled  cells 
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Figure  1 .  Velocity  models  from  straight  and  curved  ray  inversion  methods.  The  true  model  is  on  the  left.  The  inverted 
model  from  the  straight  ray  inversion  method  is  in  the  center  and  the  curved  ray  inversion  model  is  on  the  right.  The 
grayscale  indicates  velocities  in  m/ms.  The  stars  are  sources  and  the  circles  are  receivers.  White  spots  in  the  figures 
above  2  m  are  where  velocities  are  greater  than  150  m/ms,  or  where  rays  do  not  sample  the  grid,  as  along  the  bottom 
of  the  straight  ray  model. 

because  they  affect  less  of  the  travel  times.  Thus,  the  edges  of  the  model  are  less  constrained  than 
the  middle  and,  consequently,  our  uncertainty  is  greater  for  the  edges.  From  the  ray  coverage,  we 
have  more  confidence  in  the  modeled  features  in  the  region  between  3  to  7  m  depth  and  1.5  to  4  m 
distance. 

The  curved  ray  inversion  has  a  much  different  ray  coverage  character.  Thin  zones  of 
higher  coverage  exist  in  the  model.  High  ray  density  zones  extend  horizontally,  most  notably  at 
1.5  to  2.0  m.  Another  strong  horizontal  band  is  between  8.5  to  9  m.  These  bands  correspond  to  the 
high  velocity  zones  in  the  input  model.  Diagonal  zones  of  high  ray  density  also  exist  in  the  ray 
coverage.  As  in  the  straight  ray  case,  the  edges  of  the  model  are  sampled  by  fewer  rays.  Again, 
our  confidence  in  the  middle  portion  of  the  model  is  higher. 

Noting  the  ray  coverage,  we  would  still  have  confidence  in  the  models’  layering,  but  we 
would  be  less  certain  of  the  lateral  velocity  changes,  especially  at  the  model  boundaries.  We  could 
improve  the  ray  coverage  at  the  edges  by  increasing  our  cell  size,  thus  including  more  data  per 
cell.  However,  we  must  have  a  sufficiently  fine  grid  spacing  to  accurately  trace  the  rays  in  the 
finite-difference  eikonal  approximation  and  to  adequately  resolve  model  features.  Another 
method  to  improve  confidence  along  boundaries  is  to  incorporate  information  from  vertical  radar 
profiles  (Knoll  and  Clement,  1999)  and  geophysical  logging. 

The  different  ray  tracing  routines  in  the  inversions  account  for  the  different  inverted  veloc¬ 
ity  models.  Boundaries  are  sharper  and  velocity  changes  are  more  distinct  in  the  curved  ray 
model.  The  ray  coverage  in  the  curved  ray  method  shows  that  more  rays  are  concentrated  along 
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Figure  2.  The  ray  coverage  of  the  straight  ray  and  the  curved  ray  inversion.  The  grayscale  indicates  the  changes  in  ray 
density.  White  areas  are  where  rays  do  not  sufficiently  sample  the  model. 


the  boundaries  in  the  model,  especially  along  large  velocity  contrasts  and  high  velocity  zones.  The 
curved  rays  more  accurately  model  the  refraction  at  the  velocity  boundaries  since,  by  definition, 
the  straight  rays  can  not  bend. 


Constrained  Inversion:  Fixed  Velocity  Values 

A  variety  of  different  types  of  constraints  can  be  used  in  our  effort  to  improve  the  resolu¬ 
tion  and  accuracy  of  tomographic  inversions.  One  type  of  constraint  is  to  force  specific  velocity 
values  on  the  inverse  solution,  such  as  specific  layer  properties.  Tweeton  (1988)  has  shown  the 
value  of  including  layer  constraints  in  the  inverse  procedure.  We  ran  a  straight  ray  simulation, 
shown  in  Figure  3,  to  demonstrate  this  point.  Constraining  the  upper  2  meters  to  be  laterally  con¬ 
stant  produces  a  model  that  more  sharply  defines  the  water  table  boundary.  The  transition  zone 
from  the  high  velocity  upper  layer  to  the  slower  layer  beneath  is  more  sharply  defined  both  verti¬ 
cally  and  horizontally.  Moreover,  the  layering  beneath  the  water  table  is  more  clearly  defined  as 
are  the  velocity  changes  in  the  two  heterogeneous  velocity  layers.  This  relatively  simple  assump¬ 
tion,  constant  lateral  velocity  in  the  vadose  zone,  greatly  improves  the  model  reconstruction. 

Constrained  Inversion:  Regularization 

We  also  investigated  the  use  of  different  types  of  regularization  as  constraints  for  the 
curved  ray  inverse  procedure.  We  inverted  the  noise-free  data  using  a  flatness  (first  difference) 
constraint  and  a  smoothness  (second  difference)  constraint  (Aldridge  and  Oldenburg,  1993).  The 
flatness  constraint  seeks  to  find  the  model  that  does  not  contain  velocity  changes.  In  other  words, 
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Figure  3.  Effects  of  constraining  upper  layers  to  imitate  the  vadose  zone.  The  center  figure  is  inverted  with  smoothing 
only.  In  the  right  figure,  the  upper  2  m  are  constrained  to  have  a  laterally  constant  velocity.  The  grayscale  indicates  the 
velocity  in  m/jis.  Note  that  the  layering  in  the  right  figure  is  sharper  and  that  the  lateral  velocity  changes  are  more 
apparent  and  their  velocity  values  contrast  more. 

we  want  a  solution  with  as  little  change  as  possible.  The  smoothness  constraint  seeks  to  minimize 
the  gradient  of  the  velocity  change  in  the  designated  direction,  that  is,  we  want  a  solution  that  has 
small  changes  in  velocity  gradient. 

Figure  4  shows  the  results  of  the  two  different  constraints  on  the  curved  ray  inversion.  The 
velocity  variation  is  greater  horizontally  in  the  smoothness  model  than  in  the  flatness  model.  This 
variation  reflects  the  nature  of  the  constraint;  the  smoothness  constraint  limits  velocity  gradient 
changes  whereas  the  flatness  constraint  limits  velocity  changes.  The  ray  coverage  is  similar 
between  the  two  different  models.  Again,  the  similarity  between  the  models  indicates  that  the 
inversion  is  adequately  modeling  the  data  and  that  the  choice  of  regularization  constraint  is  not 
critical,  at  least  for  this  earth  model. 


Effect  of  Grid  Size 

One  decision  that  must  be  made  in  parameterizing  a  tomography  problem  is  the  size  of  the 
modeling  grid  cells.  We  used  the  straight  ray  algorithm  to  explore  this  issue.  First,  we  inverted  the 
data  on  a  0.5  m  grid.  In  figure  5,  we  compare  this  result  to  the  model  using  a  0.1  m  cell  size.  The 
0.5  m  model  contains  larger  amplitude  velocity  variations  that  are  more  sharply  defined.  Although 
the  larger  grid  spacing  better  replicates  the  velocity  variations,  the  model  contains  unreasonable 
velocity  fluctuations,  especially  at  the  edges. 

Importance  of  Starting  Model 

Another  important  parameter  for  curved  ray  tomographic  inversions  is  the  starting  model 
for  the  nonlinear  inversion.  The  inversion  process  traces  rays  through  the  starting  model,  then 
alters  the  model  to  reduce  the  misfit  error.  A  danger  in  this  iterative  method  is  that  the  routine 
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Figure  4.  The  tomogram  in  the  center  was  computed  with  a  flatness  constraint.  The  tomogram  on  the  right  was 
computed  with  a  smoothness  constraint.  The  grayscale  indicates  the  velocity  in  m/ps.  The  flatness  constraint  is  the 
first  derivative  and  the  smoothness  constraint  is  the  second  derivative.  White  regions  indicate  velocities  greater  than 
0.150  m/ps. 
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Figure  5.  Comparing  the  effects  of  grid  size  in  the  model.  The  figure  in  the  center  was  computed  using  0.1  m  grid 
spacing.  The  right  figure  has  a  grid  spacing  of  0.5  m.  The  grayscale  indicates  the  velocity  in  m/ps.  The  white  regions 
in  the  upper  2  m  represent  velocities  greater  than  0.150  m/ps.  The  white  zone  at  the  bottom  of  the  center  figure  indi¬ 
cates  a  lack  of  ray  coverage. 
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finds  a  local  minimum,  not  the  global  minimum  solution.  To  test  that  the  solution  is  near  the  glo¬ 
bal  minimum  of  the  objective  function,  we  inverted  the  data  using  several  different  starting  mod¬ 
els  consisting  of  a  homogeneous  medium  with  a  single  velocity.  We  used  velocities  of  0.05,  0.09 
and  0.20  m/ns  to  cover  the  range  of  realistic  velocities  for  ground  penetrating  radar  studies.  The 
resulting  tomograms  are  essentially  the  same,  suggesting  that  the  solution  is  close  to  the  global 
minimum. 


Effect  of  Noise 

To  investigate  the  effect  of  noise,  we  generated  and  then  inverted  two  data  sets,  the  first 
with  no  noise  and  the  second  with  noise.  We  also  added  random  noise  to  each  travel  time  to  simu¬ 
late  noisy  data.  The  noise  ranged  from  -0.5  to  0.5  ns.  The  results  for  the  curved  ray  inversion  are 
shown  in  Figure  6.  Both  methods  model  the  noisy  data  about  as  well  as  they  did  the  noise  free 
data.  The  inversion  statistics  are  very  similar  to  the  noise-free  statistics,  except  that  the  RMS  error 
values  are  larger,  not  an  unexpected  result.  The  level  of  noise  that  we  added  is  realistic  for  bore¬ 
hole  radar  data  (Saintenoy  and  Scales,  personal  communication). 


Model  Noise  free  data  Noisy  data 


Distance  (m)  Distance  (m)  Distance  (m) 


Figure  6.  Comparison  of  inverting  noise-free  and  noisy  data.  The  two  models  were  inverted  with  the  curved  ray  algo¬ 
rithm  with  identical  parameters.  The  grayscale  corresponds  to  velocities  in  m/ps.  Similar  features  are  reconstructed  in 
both  figures.  The  noisy  data  has  a  larger  variance  in  error  misfit. 


Effect  of  Ray  Coverage 

Finally,  we  study  the  effects  of  limited  ray  coverage  on  the  model.  The  full  data  set  has 
rays  that  encompass  angles  between  -70°  and  70°.  In  figure  7,  we  compare  the  full  data  set  with 
data  whose  angular  coverage  is  windowed  to  those  rays  between  -40°  and  40°  and  another  set 

between  -10°  and  10°.  The  wider  angular  coverage  better  resolves  the  lateral  horizontal  velocity 
changes  as  well  as  the  layer  boundaries.  As  the  angular  coverage  narrows,  the  lateral  resolution 

decreases.  In  the  -10°  to  10°  case,  the  model  consists  of  horizontal  layers  with  no  apparent  lateral 
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changes.  With  near  surface  tomography  acquisition,  wide  angular  coverages  are  limited  by  the 
distance  the  energy  can  propagate,  not  by  the  geometry  of  the  experiment. 


Figure  7.  Effects  of  limited  ray  coverage.  The  full  data  set  is  used  to  compute  the  tomogram  on  the  left.  The  angular 
coverage  is  between  -70°  to  70°  from  the  horizontal.  The  tomogram  in  the  center  is  computed  with  a  limits  aperture 

window  with  rays  traveling  between  -40°  and  40°,  and  the  right  tomogram  has  a  window  of  -10°  to  10°.  The  gray¬ 
scale  indicates  the  velocity  in  m/(is. 


SUMMARY 

We  have  presented  a  synthetic  modeling  study  to  examine  some  of  the  issues  that  effect 
tomographic  inversion  of  crosshole  radar  data.  We  looked  at  the  effects  of  the  forward  model,  the 
type  of  model  constraints,  the  grid  size,  and  the  starting  model.  These  factors  are  controlled  by  the 
inversion  routines.  We  also  examined  the  effect  of  noise  in  the  data  and  ray  coverage  through  the 
model.  These  effects  are  controlled  by  the  field  acquisition,  and  once  the  data  has  been  collected, 
we  can  not  change  these  factors. 

The  greatest  effect  on  the  inversion  results  is  the  choice  of  the  forward  model.  The  curved 
ray  method  better  imaged  the  vertical  and  horizontal  velocity  changes.  Another  important  effect 
on  the  inversion  is  the  structure  imposed  on  the  solution.  Inversions  with  some  prior  information 
about  the  structure  of  the  sampled  medium  were  better  able  to  reconstruct  the  heterogeneities  in 
the  medium.  Both  curved  and  straight  ray  methods  benefited  from  this  type  of  information. 
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